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In the enclosure, parts of the paper that do not relate to the invention have not 
been reproduced. Section 3.1 of the paper, which carries the heading "Aerosol 
Retrieval Method" describes the conception and reduction to practice of the 
invention of claim 1 of the subject patent application. See also equation (2) on a 
previous page (the pages are not numbered). 

The paper describes our improved method of correcting for atmospheric effects on 
a remote image of the Earth's surface taken from above, wherein the image 
comprises a number of images of the same scene each including a large number of 
pixels, each at a different wavelength band, and including infrared through visible 
wavelengths, the method comprising providing a radiation transport model that 
relates spectral radiance to spectral reflectance via a set of parameters; 
providing a discrete number of trial aerosol visibility values for at least one of one 
or both of trial aerosol property values and aerosol types; using the radiation 
transport model to calculate the model parameter values for each of the trial 
aerosol visibility values; selecting image pixels having a one or more presumed, 
predefined ratios of reflectances among two or more specific wavelength bands; 
using the radiation transport model parameters to determine the surface 
reflectance for the selected image pixels for each of the specific wavelength bands 
for each combination of trial visibility value and trial aerosol property value or 
values, or aerosol type; comparing the determined surface reflectances to the 
predefined ratio of reflectances; and resolving from the comparison a corrected 
image visibility value for each trial aerosol property value or values or aerosol 
type. 
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6. The last paragraph of Section 3.1 of the paper describes our use of the claimed 
method on data from two different imaging sensors. The reported results of our 
use of the claimed method evidences our reduction to practice before July 14, 
2000. 

7. The enclosure establishes that the inventors had conceived of and reduced to 
practice the invention of claim 1 before July 14, 2000. 
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STATUS OF ATMOSPHERIC CORRECTION USING A MODTRAN4-BASED ALGORITHM 



Michael W. Matthew, 1 Steven M. Adler-Golden, 1 Alexander Berk, 1 Steven C. Richtsmeier, 1 Robert Y. Levine, 1 
Lawrence S. Bernstein, 1 Prabhat K. Acharya, 1 Gail P. Anderson, 2 Gerry W. Felde, 2 Michael P. Hoke, 2 
Anthony Ratkowski, 2 Hsiao-Hua Burke, 3 Robert D. Kaiser, 4 and David P. Miller 4 

1. INTRODUCTION 



2. METHODOLOGY 

2.1 Review of Basic Method 

A brief review of the AFRL/SSI atmospheric correction method is presented. We start from a standard 
equation for spectral radiance at a sensor pixel, L*, that applies to the solar wavelength range (thermal emission is 
neglected) and flat, Lambertian materials or their equivalents. The equation can be written as (Vermote et ah, 1994) 

L* = Ap/( 1 - p e S) + Bp e /( 1 -p e S) + L* a ( 1 ) 

Here p is the pixel surface reflectance, p e is an average surface reflectance for the pixel and a surrounding region, S 
is the spherical albedo of the atmosphere, L* a is the radiance backscattered by the atmosphere, and A and B are 
coefficients that depend on atmospheric and geometric conditions but not on the surface. Each of these variables 
depends on the spectral channel; the wavelength index has been omitted for simplicity. The first term in Equation 
(1) corresponds to radiance that is reflected from the surface and travels directly into the sensor, while the second 
term corresponds to radiance from the surface that is scattered by the atmosphere into the sensor. The distinction 
between p and p € accounts for the "adjacency effect" (spatial mixing of radiance among nearby pixels) caused by 
atmospheric scattering. The adjacency effect correction may be ignored by setting p e = p. However, this can result 
in significant reflectance errors at short wavelengths, especially under hazy conditions and when there are strong 
contrasts among the materials in the scene (Adler-Golden et al., 1999). 
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The values of A, B, S and L* a are determined from MODTRAN4 calculations that use the viewing and solar 
angles and the mean surface elevation of the measurement and assume a certain model atmosphere, aerosol type, and 
visible range. The values of A, B, S and L* a are strongly dependent on the water vapor column amount, which is 
generally not well known and may vary across the scene. To account for unknown and variable column water vapor, 
the MODTRAN4 calculations are looped over a series of different column amounts, then selected wavelength 
channels of the image are analyzed to retrieve an estimated amount for each pixel. Specifically, radiance averages 
are gathered for two sets of channels, an "absorption" set centered at a water band (typically the 1.13 urn band) and a 
"reference" set of channels taken from just outside the band. A 2-dimensional look-up table (LUT) for retrieving the 
water vapor from these radiances is constructed. One dimension of the table is the reference to absorption ratio and 
the other is the reference radiance. The second dimension accounts for a reflectance-dependent variation in the ratio 
arising from the different amounts of absorption in the atmospherically-scattered and surface-reflected components 
of the radiance. After the water retrieval is performed, Equation (1) is solved for the pixel surface reflectances in all 
of the sensor channels. The solution method (Richter, 1996; Vermote et ai, 1997) involves computing a spatially 
averaged radiance image Z,* e , from which the spatially averaged reflectance p e is estimated using the approximate 
equation 

L* e ^(A+B) P( /(l-p e S)+L\ (2) 

The spatial averaging is performed using a point-spread function that describes the relative contributions to the pixel 
radiance from points on the ground at different distances from the direct line of sight. The SSI/AFRL algorithm 
approximates this function as a radial exponential. 

In the above discussion it has been assumed that the quantity of aerosol or haze in the scene has been 
adequately estimated. As described previously (Adler-Golden et a!., 1999), the SSI/AFRL algorithm includes a 
method for retrieving an estimated aerosol/haze amount from one or more reference surfaces in the scene that have a 
known reflectance in some wavelength bandpass. Best results are obtained using short (visible) wavelengths and 
either a very dark surface, such as vegetation or deep calm water, or a very bright surface, such as a white calibration 
target that is large enough to fill a whole pixel. In this method, calculations to determine A, B, S and L* a are carried 
out for the spectral channels in the designated bandpass. Instead of iterating over different water vapor values, these 
calculations are performed over a series of visible ranges, e.g. 200, 100, 50, 33, 25, 20 and 17 km, that are evenly 
spaced in their reciprocals (optical depths). The user selects the reference pixels and assigns them a mean 
reflectance value for the selected channels. The algorithm derives a visible range for each reference pixel by 
interpolating from a 2-D LUT that depends on L and L* e . From these results an average or other "best" estimate of 
the visible range can be derived and used for the MODTRAN4 calculation loop over water vapor. 

An example of AVIRIS data processed with the above procedure is shown in Figure 1 . The data are of 
white and black calibration panels at the Stennis Space Center, and were acquired from the low altitude (Twin Otter) 
platform in October, 1998, After adjusting the wavelength calibration slightly and performing the atmospheric 
correction, the spectra were smoothed using a "polishing" algorithm (Adler-Golden et al, 1999; Boardman, 1998). 
The AVIRIS spectra (particularly the white panel) show some absorption residuals adjacent to the cut-out regions of 
very strong water absorption and also at the 0.94 urn water band, but on the whole the agreement with the "ground 
truth" spectra is good. The -0.01 difference between the black panel spectra is within the variability and uncertainty 
in the ground truth measurement. 

2.2 Limitations 

The basic atmospheric correction method described above, as well as those in other first-principles codes 
such as ATREM, work well in many but not all scenes. In particular, they require cloud-free conditions, the 
presence of at least one material in the scene with a known reflectance at a visible wavelength, and sufficient 
computing time to perform tens of mathematical operations per image pixel per wavelength channel. 
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Figure 1 . Comparisons of Atmospherically Corrected AVIRJS Data with Ground Truth Measurements for White 
(left) and Black (right) Calibration Panels at the Stennis Space Center. 

Cloud Effects. Clouds and cloud shadows pose several problems for atmospheric correction. Not only do 
cloud-contaminated pixels have incorrect reflectance, they also can degrade the reflectance accuracy in other parts of 
the scene. This is because clouds impact the spatially averaged radiance L* e used in Equation (2) to generate p e for 
the adjacency correction. According to theory, p e should account only for reflecting material that is below the 
scattering atmosphere. While clouds typically lie below most of the molecular (Rayleigh) scattering, which is 
important at the very shortest (blue-violet) wavelengths, they typically lie above the aerosol and haze layers that 
dominate the scattering in the rest of the spectrum. Since clouds are typically much brighter than the terrain, p e is 
overestimated, leading to underestimated surface reflectance retrievals. Therefore it is important to identify and 
remove cloud-contaminated pixels prior to the calculation of p e . Cloud shadows also produce reflectance errors; 
however, their effects on material identification can be compensated to some extent, and their impact on clear parts 
of the scene is minimal. 



Absence of Accurately Known Surfaces. For the purpose of aerosol/haze amount retrieval, vegetation, 
water, or other dark surfaces can frequently be identified in a scene. However, reflectance values for these surfaces 
at appropriate wavelengths are often not known to within the accuracy needed (around 0.01 reflectance units or 
better) for a good retrieval. Even with "calibrated" surfaces the reflectance may not be known to within this 
accuracy because of complications caused by non-Lambertian bidirectional reflectance distribution functions. 
However, a method based on a known reflectance ratio for different wavelengths, such as the (Kaufman et al. 1997) 
dark pixel method, can minimize these problems. 

Computing Time Requirements. For a typical image containing several hundred or more spectral channels 
and hundreds of thousands of pixels or more, the speed of the atmospheric correction is fundamentally limited by the 
mathematical operations required to generate the reflectance values for each pixel and channel from the Equation (1) 
parameters. Current algorithms, such as ATREM and the SSI/AFRL code, that use pixel-specific values of water 
vapor (and possibly other quantities such as p e ) require tens of operations per pixel-channel. Most of the operations 
are consumed in interpolating to find the appropriate A t B, S and L* a parameters for each pixel. A more efficient 
procedure is needed to achieve high-speed atmospheric correction. 

Model Accuracy. Any first-principles atmospheric correction method is necessarily limited by the accuracy 
of its radiation transport model. We have recently had opportunities to validate MODTRAN4 against analytical and 
Monte Carlo scattering calculations as well as against "exact" line-by-line transmittance calculations degraded to 
AVIRIS spectral resolution. Excellent agreement was obtained in each case. However, consistency among 
calculations does not guarantee an accurate representation of reality. Recently, it has been discovered that a number 
of bands of water vapor in the HITRAN atlas (Rothman et al., 1998), upon which MODTRAN4 as well as line-by- 
line codes are based, have incorrect line strengths. In the 0.94 urn band the errors are around 14%, and signficantly 
impact the atmospheric correction. These errors have been corrected in the most recent version of MODTRAN4. 



Other known deficiencies in MODTRAN4 include the omission of certain collision-induced absorption bands of 
oxygen (Solomon et al, 1998). 

3. UPGRADES 

The AFRL/SSI atmospheric correction code has been upgraded with several new algorithms that address, if 
not completely solve, the abovementioned limitations: 

• A new method has been implemented for retrieving the aerosol/haze amount from an assumed ratio of 
in-band reflectances, rather than from an assumed reflectance value. This method can utilize user- 
selected pixels or can automatically find suitable dark terrain pixels (Kaufman et al., 1997) for the 
retrieval. 

• An algorithm for identifying cloud-containing pixels in AVIRIS or similar data has been implemented, 
and is used to improve the calculation of L* e and p e in Equations (1) and (2). Since this algorithm 
requires a prior water vapor retrieval, the order and number of steps in the atmospheric correction 
process has been altered. 

• A new method has been developed that greatly reduces the number of mathematical operations required 
to generate the reflectance values once the atmospheric parameters have been defined. The method 
operates by averaging the water vapor and p e values over small groups of neighboring pixels, so that 
the same A, B, S, and L* a parameters may be assigned to all pixels in the group. 

These upgrades are described in more detail below. In addition, we show an example of the impact of the water line 
strength revisions on a reflectance retrieval. 

3.1 Aerosol Retrieval Method 

A general reflectance ratio-based algorithm has been developed for retrieving an aerosol amount (i.e., the 
visible range). The reference pixels can be chosen by the user, or dark pixels can be selected automatically based on 
a specified maximum reflectance. To implement the (Kaufman et ai 1997) method, one chooses bandpasses 
centered at 0.66 um and 2. 1 u,m, a reflectance ratio of -0.5, and a 2. 1 \im reflectance maximum of around 0. 1 . 

Radiance images in each of the two bandpasses are assembled from both the original data cube and from the 
spatially averaged radiance L* e . MODTRAN4 calculations are conducted to determine A, B, S and L* a for a series 
of trial visible range values. For each visible range and reference pixel, the reflectance solutions for the two 
bandpasses are calculated, and the reflectance error for the shorter-wavelength bandpass (the difference between the 
calculated reflectance and the calculated longer-wavelength reflectance times the assumed ratio) is tabulated. A 
visible range estimate for each selected pixel can be obtained by interpolating within the resulting array of 
reflectance errors to find the value that yields zero error. To more efficiently calculate a scene-average visible range, 
the reflectance error arrays are averaged over all reference pixels, and the interpolation is performed on the result. 

We have tested this method using different reference materials and data from two different imaging sensors, 
including AVIRIS. Using calibration panels as reference pixels, the visible range results were compared with results 
from the original reflectance-based method, and very good agreement was found. Using natural dark terrain, results 
were assessed for different reflectance cutoffs and ratio values within the tolerances found by (Kaufman et ai 1997). 
Our preliminary analysis indicates that the typical obtainable retrieval accuracy is 0.01 to 0.02 per km for Invisible 
range). For example, the difference between retrieved visible ranges values of 50 km and 300 km may not be 
significant, whereas the differences between values of 20 km and 50 km or 13 km and 20 km would be considered 
significant. In some cases the results can be made less sensitive to the value of the reflectance ratio by choosing a 
very low reflectance cutoff, such as 0.04, for the dark pixel selection. However, for scenes containing shallow or 
turbid water bodies a low cutoff can produce anomalous results, since the low cutoff favors the selection of the water 
pixels, which can have a very different reflectance ratio than dark land pixels. A more sophisticated pixel selection 



method, such as one that includes both maximum and minimum reflectance cutoffs to discriminate against surface 
water, should provide better results. 

3.2 Identification and Utilization of Cloud-Containing Pixels 

An algorithm has been developed for generating a cloud "mask" that identifies cloud-containing pixels in 
the scene. At the present time the main uses of the cloud mask in the SSI/AFRL atmospheric correction code are (1) 
to indicate regions where the atmospheric correction is invalid or suspect, and (2) to flag pixels that need to be 
removed from the calculation of L* e (currently these pixels are replaced by the scene average radiance). While it is 
most important to flag bright, opaque clouds, it is also desirable to find pixels that contain appreciable cirrus or other 
thin cloud cover. 

A number of cloud detection algorithms have been developed based on multispectral data. A 
comprehensive review is presented by Ackerman et al. (1998), who developed an algorithm for the MODIS sensor. 
Their algorithm uses a combination of tests, including (1) a color balance test based on a SWIR/red reflectance ratio 
(0.9 < p(0.87 urn) /p(0.66 um) < 1 . 1 indicates clouds], (2) a reflectance test at 0.94 um (a high signal correlates with 
low column water vapor, hence reflection from a bright, elevated object), and (3) a variety of brightness tests at IR 
wavelengths. In a paper on simple algorithms for multispectral atmospheric correction, (Borel et al 1999) discuss 
an analogue to the SWIR/red ratio test that combines an upper threshold on the NDVI (Normalized Differential 
Vegetation Index) with a lower threshold on the SWIR bandpass (i.e., a brightness test). They also describe a water 
vapor absorption test involving a continuum interpolated band ratio (CIBR). 

With AVIRJS there is no IR coverage past 2.5 um, but there is high spectral resolution that permits a very 
good column water vapor retrieval. Therefore, following Ackerman and Borel, we devised a cloud mask based on 
combining tests for brightness, color balance, and low column water in the visible and SWIR regions. Because of 
processing time constraints, it is advantageous to utilize bands that are already being gathered by the atmospheric 
correction code for other purposes (e.g., 2.1 um, 1.13 u.m water absorption and reference bandpasses, 0.66 um, red, 
blue, and green bandpasses used for image display) and to use the retrieved water vapor amounts. From these data 
we implemented analogues of the tests described above. 

From the standpoint of clear-sky atmospheric correction, the main effect of clouds in the scene arises from 
the adjacency effect compensation, which requires a spatially smoothed radiance. It is not appropriate to include 
cloudy pixels in the smoothed radiance, which means that the cloud mask must be determined prior to both the 
aerosol retrieval and the atmospheric correction. A cloud test based on water vapor must use some assumed aerosol 
amount, and a test based on reflectances cannot include the adjacency effect compensation. Given these 
requirements, the preferred sequence of steps for the atmospheric correction process is as follows: 

1 . Initial water vapor retrieval. A nominal visible range (e.g., 50 km) is assumed. 

2. Cloud mask generation. Brightness and color balance tests are applied to establish probable clear 
pixels, and a spatially average water vapor average is taken. Pixels containing significantly lower water 
vapor than this spatial average are identified, and the results of this test and the other tests are 
combined to define the opaque cloud mask. 

3. Spatial averaging of the radiance using the adjacency effect point-spread function. Prior to averaging, 
the scene-average radiance replaces the actual radiance in the cloud-masked pixels. 

4. Aerosol (visible range) retrieval. The automated ratio-based algorithm is used with adjacency 
correction (both the smoothed and unsmoothed radiances are input). 



5. Refined water vapor retrieval. The derived visible range and perhaps a narrowed range of water 
column amounts are used.O 

6. The cloud mask may be recalculated, but it should not be much different than before. 

7. Full reflectance spectrum retrieval. 

A convenient method for the cloud mask generation, incorporated in the most recent version of the 
SSI/AFRL atmospheric correction method, is outlined below. Brightness, color balance, and water vapor tests are 
used together to define a mask for "ordinary", low-altitude clouds. In addition, 1 .38 urn data are used to define a 
separate mask for high-altitude (i.e., cirrus) clouds, following the work of Gao and co-workers (1998, 1993). 

The brightness test requires that an atmospheric correction from radiance to reflectance units be performed 
for at least one sensor bandpass. Since the water reference reflectance channel average (taken from either side of 
1.13 um) and a corresponding reflectance are already generated, it can be used for the brightness test. Borel et al 
(1999) recommend a reflectance lower threshold of around 0.3 for clouds in the SWIR. We have obtained good 
results with a similar value, 0.4. 

The color balance test involves comparing at least two bandpasses at different wavelengths. One bandpass 
can be the water reference, the second can be a visible bandpass, preferably green wavelengths, properly scaled. The 
test outcome is positive if the ratio of effective reflectances (radiance divided by the solar function) in the green and 
water' reference bandpasses is unity to within some bounding values. Suitable bounding values determined by trial 
and error are 0.4 and 1 .2. 

The low-water test involves comparing the pixel's column water vapor with a threshold value that is derived 
from pixels that fail both the brightness and color balance tests and therefore are classified as clear. The threshold is 
defined with respect to a clear-pixel spatial average, obtained by convolving the clear pixel image with a window 
that is smaller than the image. Ideally, the window should be larger than typical cloud dimensions but smaller than 
typical large-scale topographic dimensions. For example, for AVIRIS data taken from a 20 km altitude, a suitable 
window size is around 40 x 40 pixels. The outcome of the low- water test is defined to be positive for a pixel if its 
column water vapor is less than 85% of the clear-pixel value. 

To generate the high cloud mask, one or two channels of data in the center of the 1 .38 um water band are 
selected. The data are histogrammed, the maximum of the histogram is assigned to the background level, and pixels 
whose signals exceed some threshold (presently 0.03 uW/nm/cm 2 /sr) above the background level are flagged. We 
have found that this method often detects thick, lower-altitude clouds as well as cirrus clouds. 

To date, limited testing of the cloud identification algorithm has been conducted. An application to an 
AVIRIS image taken near North Conway, NH is shown in Figure 2. The pixels are color coded according to which 
cloud tests came out positive. Except where the clouds are extremely thin, all pixels that appear to the eye to be 
contaminated with clouds are flagged by the algorithm, and false positives are not evident. A low rate of false 
positives, particularly for the "ordinary" cloud test, has been verified using a variety of cloud-free scenes. 



Figure 2. Left, RGB Rendition of AVIRIS Radiance Image, Right, calculated cloud mask. The pixels tested 
positive as follows: Yellow = high cloud, Violet = "ordinary" cloud, Red = high + "ordinary". Green pixels tested 
positive for high clouds and also had column water values below the lower limit of the water LUT, indicative of a 
thick cloud. 

3.3 Reflectance Calculation Speedup 

The reflectance calculation described in Section 2. 1 can be made much faster, with little sacrifice in 
accuracy, by approximating A, B, S, L* a and p e using average values for a group of nearby pixels (such as an N x N 
array), referred to here as a "superpixel". The method takes advantage of the fact that Equation (1) relating radiance 
to reflectance can be transformed into the simple linear equation, 

p = mL* + b (3) 

where m and b are expressed in terms of superpixel values for .4, B, S, L* a , and p e . Suitable values for these 
parameters are determined by using superpixel-average water vapor amounts to interpolate from MODTRAN4- 
derived LUTs. The superpixel water vapor amounts may be either averages of the retrievals from individual pixels 
or retrievals from superpixel-average radiances. p e is calculated from the value of L* e for the superpixel. Note that 
since L*. is itself a spatial average, for all practical purposes it does not need to be calculated on a single-pixel basis 
in the first place. Once m and b are defined, p is calculated from L* for each pixel. 

The speedup in the calculation of the reflectance compared to the standard pixel-by-pixel approach derives 
from the fact that the interpolations and other mathematical steps required to generate m and b (approximately 2 1 
arithmetic operations) are performed only once per N x N pixels. In the limit of large N, the number of operations 
per pixel-channel reduces to the 2 operations in Equation (3), which are the same as in the Empirical Line Method. 
Most of the speed benefit can be achieved even with a modest superpixel size, such as N = 4 (see Table 1), which 
yields only marginal differences with the "exact" N = 1 results. 

To date, we have implemented the superpixel method in an IDL language code and obtained a fourfold 
improvement in speed, to around 1/3 s per line of 614 AVIRIS pixels on a 330 MHz PC. A further order-of- 
magnitude speedup is anticipated with recoding to a more efficient language such as C or FORTRAN. 



Table 1. Theoretical Number of Floating Point Operations Per Pixel-channel using the Superpixel Method 
for Calculating Spectral Reflectance. 



Superpixel 
Dimensions 


FLOPs / pixel / 
channel 


lxl 


23 


2x2 
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3.4 Revised Water Line Strengths 

A revised HITRAN line list (Rothman et al, 1998), containing the (Giver et ai 2000) water line parameter 
corrections, was used to formulate a new set of water band model parameters for MODTRAN4. The impact of the 
new parameters on AVIRIS data is shown in Figure 3 for the Stennis Space Center white panel. To best show the 
residual errors, no spectral "polishing" was applied. The new parameters virtually eliminate the anomalous 0,94 um 
absorption which was found earlier, and which has been a persistent artifact in AVIRIS retrievals, especially in moist 
atmospheres. At most other wavelengths the two sets of results are nearly identical. The new parameters led to a 
very small increase in the retrieved water column amount, from 1550 to 1570 atm-cm, due to a ~1% change in the 
1 .13 um band strength. Both results are in remarkable (perhaps partly fortuitous) agreement with the value of 1560 
atm-cm measured by a radiosonde near the time and location of the AVIRIS flight. 
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Figure 3. Comparison of White Panel Reflectance Spectra retrieved from AVIRIS Data Using MODTRAN4 with 
Corrected (blue) and Original (red) HITRAN Water Line Parameters. 
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